#include "PrincipalAxisOptimizerProvider.h"

struct
{
  doublereal fx, ldt, dmin__;
  integer nf, nl;
} global_;

#define global_1 global_

struct
{
  doublereal v[400] /* was [20][20] */, q0[20], q1[20], qa, qb, qc, qd0, qd1,
    qf1;
} q_;

#define q_1 q_

/* Table of constant values */

static integer c__1 = 1;
static integer c__2 = 2;
static logical c_false = FALSE_;
static integer c__10 = 10;
static integer c__4 = 4;
static logical c_true = TRUE_;
static integer c__3 = 3;
static integer c__0 = 0;

doublereal SmartMathLibrary::InternalRoutines::PrincipalAxisOptimizerProvider
  ::PrincipalAxisMethod(doublereal *t0, doublereal *machep, doublereal *h0,
  integer *n, integer *prin, doublereal *x, PrincipalAxisHelperFunction ^ f,
  doublereal *fmin)
{
  /* System generated locals */
  integer i__1, i__2, i__3;
  doublereal ret_val, d__1;

  /* Local variables */
  static doublereal d__[20], h__;
  static integer i__, j, k;
  static doublereal s, t, y[20], z__[20], f1;
  static integer k2;
  static doublereal m2, m4, t2, df, dn;
  static integer kl, ii;
  static doublereal sf;
  static integer kt;
  static doublereal sl;
  static integer im1, km1;
  static doublereal dni, lds;
  static integer ktm;
  static doublereal scbd;
  static integer idim;
  static logical illc;
  static integer klmk;
  static doublereal ldfac, large, small, value;
  static doublereal vlarge;
  static doublereal vsmall;

  /*                             LAST MODIFIED 3/1/73 */

  /*     PRAXIS RETURNS THE MINIMUM OF THE FUNCTION F(X,N) OF N VARIABLES */
  /*     USING THE PRINCIPAL AXIS METHOD.  THE GRADIENT OF THE FUNCTION IS */
  /*     NOT REQUIRED. */

  /*     FOR A DESCRIPTION OF THE ALGORITHM, SEE CHAPTER SEVEN OF */
  /*     "ALGORITHMS FOR FINDING ZEROS AND EXTREMA OF FUNCTIONS WITHOUT */
  /*     CALCULATING DERIVATIVES" BY RICHARD P BRENT. */

  /*     THE PARAMETERS ARE: */
  /*     T0       IS A TOLERANCE.  PRAXIS ATTEMPTS TO RETURN PRAXIS=F(X) */
  /*              SUCH THAT IF X0 IS THE TRUE LOCAL MINIMUM NEAR X, THEN */
  /*              NORM(X-X0) < T0 + SQUAREROOT(MACHEP)*NORM(X). */
  /*     MACHEP   IS THE MACHINE PRECISION, THE SMALLEST NUMBER SUCH THAT */
  /*              1 + MACHEP > 1.  MACHEP SHOULD BE 16.**-13 (ABOUT */
  /*              2.22D-16) FOR REAL*8 ARITHMETIC ON THE IBM 360. */
  /*     H0       IS THE MAXIMUM STEP SIZE.  H0 SHOULD BE SET TO ABOUT THE */
  /*              MAXIMUM DISTANCE FROM THE INITIAL GUESS TO THE MINIMUM. */
  /*              (IF H0 IS SET TOO LARGE OR TOO SMALL, THE INITIAL RATE OF */
  /*              CONVERGENCE MAY BE SLOW.) */
  /*     N        (AT LEAST TWO) IS THE NUMBER OF VARIABLES UPON WHICH */
  /*              THE FUNCTION DEPENDS. */
  /*     PRIN     CONTROLS THE PRINTING OF INTERMEDIATE RESULTS. */
  /*              IF PRIN=0, NOTHING IS PRINTED. */
  /*              IF PRIN=1, F IS PRINTED AFTER EVERY N+1 OR N+2 LINEAR */
  /*              MINIMIZATIONS.  FINAL X IS PRINTED, BUT INTERMEDIATE X IS */
  /*              PRINTED ONLY IF N IS AT MOST 4. */
  /*              IF PRIN=2, THE SCALE FACTORS AND THE PRINCIPAL VALUES OF */
  /*              THE APPROXIMATING QUADRATIC FORM ARE ALSO PRINTED. */
  /*              IF PRIN=3, X IS ALSO PRINTED AFTER EVERY FEW LINEAR */
  /*              MINIMIZATIONS. */
  /*              IF PRIN=4, THE PRINCIPAL VECTORS OF THE APPROXIMATING */
  /*              QUADRATIC FORM ARE ALSO PRINTED. */
  /*     X        IS AN ARRAY CONTAINING ON ENTRY A GUESS OF THE POINT OF */
  /*              MINIMUM, ON RETURN THE ESTIMATED POINT OF MINIMUM. */
  /*     F(X,N)   IS THE FUNCTION TO BE MINIMIZED.  F SHOULD BE A REAL*8 */
  /*              FUNCTION DECLARED EXTERNAL IN THE CALLING PROGRAM. */
  /*     FMIN     IS AN ESTIMATE OF THE MINIMUM, USED ONLY IN PRINTING */
  /*              INTERMEDIATE RESULTS. */
  /*     THE APPROXIMATING QUADRATIC FORM IS */
  /*              Q(X') = F(X,N) + (1/2) * (X'-X)-TRANSPOSE * A * (X'-X) */
  /*     WHERE X IS THE BEST ESTIMATE OF THE MINIMUM AND A IS */
  /*              INVERSE(V-TRANSPOSE) * D * INVERSE(V) */
  /*     (V(*,*) IS THE MATRIX OF SEARCH DIRECTIONS; D(*) IS THE ARRAY */
  /*     OF SECOND DIFFERENCES).  IF F HAS CONTINUOUS SECOND DERIVATIVES */
  /*     NEAR X0, A WILL TEND TO THE HESSIAN OF F AT X0 AS X APPROACHES X0. */

  /*     IT IS ASSUMED THAT ON FLOATING-POINT UNDERFLOW THE RESULT IS SET */
  /*     TO ZERO. */
  /*     THE USER SHOULD OBSERVE THE COMMENT ON HEURISTIC NUMBERS AFTER */
  /*     THE INITIALIZATION OF MACHINE DEPENDENT NUMBERS. */


  /* .....IF N>20 OR IF N<20 AND YOU NEED MORE SPACE, CHANGE '20' TO THE */
  /*     LARGEST VALUE OF N IN THE NEXT CARD, IN THE CARD 'IDIM=20', AND */
  /*     IN THE DIMENSION STATEMENTS IN SUBROUTINES MINFIT,MIN,FLIN,QUAD. */


  /* .....INITIALIZATION..... */
  /*     MACHINE DEPENDENT NUMBERS: */

  /* Parameter adjustments */
  --x;

  /* Function Body */
  small =  *machep **machep;
  vsmall = small * small;
  large = 1. / small;
  vlarge = 1. / vsmall;
  m2 = System::Math::Sqrt(*machep);
  m4 = System::Math::Sqrt(m2);

  /*     HEURISTIC NUMBERS: */
  /*     IF THE AXES MAY BE BADLY SCALED (WHICH IS TO BE AVOIDED IF */
  /*     POSSIBLE), THEN SET SCBD=10.  OTHERWISE SET SCBD=1. */
  /*     IF THE PROBLEM IS KNOWN TO BE ILL-CONDITIONED, SET ILLC=TRUE. */
  /*     OTHERWISE SET ILLC=FALSE. */
  /*     KTM IS THE NUMBER OF ITERATIONS WITHOUT IMPROVEMENT BEFORE THE */
  /*     ALGORITHM TERMINATES.  KTM=4 IS VERY CAUTIOUS; USUALLY KTM=1 */
  /*     IS SATISFACTORY. */

  scbd = 1.;
  illc = FALSE_;
  ktm = 1;

  ldfac = .01;
  if (illc)
  {
    ldfac = .1;
  }
  kt = 0;
  global_1.nl = 0;
  global_1.nf = 1;
  global_1.fx = f(&x[1], n);
  q_1.qf1 = global_1.fx;
  t = small + abs(*t0);
  t2 = t;
  global_1.dmin__ = small;
  h__ =  *h0;
  if (h__ < t *100)
  {
    h__ = t * 100;
  }
  global_1.ldt = h__;
  /* .....THE FIRST SET OF SEARCH DIRECTIONS V IS THE IDENTITY MATRIX..... */
  i__1 =  *n;
  for (i__ = 1; i__ <= i__1; ++i__)
  {
    i__2 =  *n;
    for (j = 1; j <= i__2; ++j)
    {
      /* L10: */
      q_1.v[i__ + j * 20-21] = 0.;
    }
    /* L20: */
    q_1.v[i__ + i__ * 20-21] = 1.;
  }
  d__[0] = 0.;
  q_1.qd0 = 0.;
  i__1 =  *n;
  for (i__ = 1; i__ <= i__1; ++i__)
  {
    q_1.q0[i__ - 1] = x[i__];
    /* L30: */
    q_1.q1[i__ - 1] = x[i__];
  }
  /*if (*prin > 0) {
  print_(n, &x[1], prin, fmin);
  }*/

  /* .....THE MAIN LOOP STARTS HERE..... */
  L40: sf = d__[0];
  d__[0] = 0.;
  s = 0.;

  /* .....MINIMIZE ALONG THE FIRST DIRECTION V(*,1). */
  /*     FX MUST BE PASSED TO MIN BY VALUE. */
  value = global_1.fx;
  min_(n, &c__1, &c__2, d__, &s, &value, &c_false, f, &x[1], &t, machep, &h__);
  if (s > 0.)
  {
    goto L50;
  }
  i__1 =  *n;
  for (i__ = 1; i__ <= i__1; ++i__)
  {
    /* L45: */
    q_1.v[i__ - 1] =  - q_1.v[i__ - 1];
  }
  L50: if (sf > d__[0] *.9 && sf *.9 < d__[0])
  {
    goto L70;
  }
  i__1 =  *n;
  for (i__ = 2; i__ <= i__1; ++i__)
  {
    /* L60: */
    d__[i__ - 1] = 0.;
  }

  /* .....THE INNER LOOP STARTS HERE..... */
  L70: i__1 =  *n;
  for (k = 2; k <= i__1; ++k)
  {
    i__2 =  *n;
    for (i__ = 1; i__ <= i__2; ++i__)
    {
      /* L75: */
      y[i__ - 1] = x[i__];
    }
    sf = global_1.fx;
    if (kt > 0)
    {
      illc = TRUE_;
    }
    L80: kl = k;
    df = 0.;

    /* .....A RANDOM STEP FOLLOWS (TO AVOID RESOLUTION VALLEYS). */
    /*     PRAXIS ASSUMES THAT RANDOM RETURNS A RANDOM NUMBER UNIFORMLY */
    /*     DISTRIBUTED IN (0,1). */

    if (!illc)
    {
      goto L95;
    }
    i__2 =  *n;
    for (i__ = 1; i__ <= i__2; ++i__)
    {
      s = (global_1.ldt *.1 + t2 * pow_ii(&c__10, &kt))*(random_(n) - .5);
      z__[i__ - 1] = s;
      i__3 =  *n;
      for (j = 1; j <= i__3; ++j)
      {
        /* L85: */
        x[j] += s * q_1.v[j + i__ * 20-21];
      }
      /* L90: */
    }
    global_1.fx = f(&x[1], n);
    ++global_1.nf;

    /* .....MINIMIZE ALONG THE "NON-CONJUGATE" DIRECTIONS V(*,K),...,V(*,N) */

    L95: i__2 =  *n;
    for (k2 = k; k2 <= i__2; ++k2)
    {
      sl = global_1.fx;
      s = 0.;
      value = global_1.fx;
      min_(n, &k2, &c__2, &d__[k2 - 1], &s, &value, &c_false, f, &x[1], &t,
        machep, &h__);
      if (illc)
      {
        goto L97;
      }
      s = sl - global_1.fx;
      goto L99;
      L97: 
      /* Computing 2nd power */
      d__1 = s + z__[k2 - 1];
      s = d__[k2 - 1]*(d__1 *d__1);
      L99: if (df > s)
      {
        goto L105;
      }
      df = s;
      kl = k2;
      L105: ;
    }
    if (illc || df >= (d__1 =  *machep * 100 * global_1.fx, abs(d__1)))
    {
      goto L110;
    }

    /* .....IF THERE WAS NOT MUCH IMPROVEMENT ON THE FIRST TRY, SET */
    /*     ILLC=TRUE AND START THE INNER LOOP AGAIN..... */

    illc = TRUE_;
    goto L80;
    L110: 
    /*if (k == 2 && *prin > 1) {
    vcprnt_(&c__1, d__, n);
    }*/

    /* .....MINIMIZE ALONG THE "CONJUGATE" DIRECTIONS V(*,1),...,V(*,K-1) */

    km1 = k - 1;
    i__2 = km1;
    for (k2 = 1; k2 <= i__2; ++k2)
    {
      s = 0.;
      value = global_1.fx;
      min_(n, &k2, &c__2, &d__[k2 - 1], &s, &value, &c_false, f, &x[1], &t,
        machep, &h__);
      /* L120: */
    }
    f1 = global_1.fx;
    global_1.fx = sf;
    lds = 0.;
    i__2 =  *n;
    for (i__ = 1; i__ <= i__2; ++i__)
    {
      sl = x[i__];
      x[i__] = y[i__ - 1];
      sl -= y[i__ - 1];
      y[i__ - 1] = sl;
      /* L130: */
      lds += sl * sl;
    }
    lds = System::Math::Sqrt(lds);
    if (lds <= small)
    {
      goto L160;
    }

    /* .....DISCARD DIRECTION V(*,KL). */
    /*     IF NO RANDOM STEP WAS TAKEN, V(*,KL) IS THE "NON-CONJUGATE" */
    /*     DIRECTION ALONG WHICH THE GREATEST IMPROVEMENT WAS MADE..... */

    klmk = kl - k;
    if (klmk < 1)
    {
      goto L141;
    }
    i__2 = klmk;
    for (ii = 1; ii <= i__2; ++ii)
    {
      i__ = kl - ii;
      i__3 =  *n;
      for (j = 1; j <= i__3; ++j)
      {
        /* L135: */
        q_1.v[j + (i__ + 1) *20-21] = q_1.v[j + i__ * 20-21];
      }
      /* L140: */
      d__[i__] = d__[i__ - 1];
    }
    L141: d__[k - 1] = 0.;
    i__2 =  *n;
    for (i__ = 1; i__ <= i__2; ++i__)
    {
      /* L145: */
      q_1.v[i__ + k * 20-21] = y[i__ - 1] / lds;
    }

    /* .....MINIMIZE ALONG THE NEW "CONJUGATE" DIRECTION V(*,K), WHICH IS */
    /*     THE NORMALIZED VECTOR:  (NEW X) - (0LD X)..... */

    value = f1;
    min_(n, &k, &c__4, &d__[k - 1], &lds, &value, &c_true, f, &x[1],  &t,
      machep, &h__);
    if (lds > 0.)
    {
      goto L160;
    }
    lds =  - lds;
    i__2 =  *n;
    for (i__ = 1; i__ <= i__2; ++i__)
    {
      /* L150: */
      q_1.v[i__ + k * 20-21] =  - q_1.v[i__ + k * 20-21];
    }
    L160: global_1.ldt = ldfac * global_1.ldt;
    if (global_1.ldt < lds)
    {
      global_1.ldt = lds;
    }
    /*if (*prin > 0) {
    print_(n, &x[1], prin, fmin);
    }*/
    t2 = 0.;
    i__2 =  *n;
    for (i__ = 1; i__ <= i__2; ++i__)
    {
      /* L165: */
      /* Computing 2nd power */
      d__1 = x[i__];
      t2 += d__1 * d__1;
    }
    t2 = m2 * System::Math::Sqrt(t2) + t;

    /* .....SEE WHETHER THE LENGTH OF THE STEP TAKEN SINCE STARTING THE */
    /*     INNER LOOP EXCEEDS HALF THE TOLERANCE..... */

    if (global_1.ldt > t2 *(float).5)
    {
      kt =  - 1;
    }
    ++kt;
    if (kt > ktm)
    {
      goto L400;
    }
    /* L170: */
  }
  /* .....THE INNER LOOP ENDS HERE. */

  /*     TRY QUADRATIC EXTRAPOLATION IN CASE WE ARE IN A CURVED VALLEY. */

  /* L171: */
  quad_(n, f, &x[1], &t, machep, &h__);
  dn = 0.;
  i__1 =  *n;
  for (i__ = 1; i__ <= i__1; ++i__)
  {
    d__[i__ - 1] = 1. / System::Math::Sqrt(d__[i__ - 1]);
    if (dn < d__[i__ - 1])
    {
      dn = d__[i__ - 1];
    }
    /* L175: */
  }
  /* if (*prin > 3) {
  maprnt_(&c__1, q_1.v, &idim, n);
  }*/
  i__1 =  *n;
  for (j = 1; j <= i__1; ++j)
  {
    s = d__[j - 1] / dn;
    i__2 =  *n;
    for (i__ = 1; i__ <= i__2; ++i__)
    {
      /* L180: */
      q_1.v[i__ + j * 20-21] = s * q_1.v[i__ + j * 20-21];
    }
  }

  /* .....SCALE THE AXES TO TRY TO REDUCE THE CONDITION NUMBER..... */

  if (scbd <= 1.)
  {
    goto L200;
  }
  s = vlarge;
  i__2 =  *n;
  for (i__ = 1; i__ <= i__2; ++i__)
  {
    sl = 0.;
    i__1 =  *n;
    for (j = 1; j <= i__1; ++j)
    {
      /* L182: */
      sl += q_1.v[i__ + j * 20-21] *q_1.v[i__ + j * 20-21];
    }
    z__[i__ - 1] = System::Math::Sqrt(sl);
    if (z__[i__ - 1] < m4)
    {
      z__[i__ - 1] = m4;
    }
    if (s > z__[i__ - 1])
    {
      s = z__[i__ - 1];
    }
    /* L185: */
  }
  i__2 =  *n;
  for (i__ = 1; i__ <= i__2; ++i__)
  {
    sl = s / z__[i__ - 1];
    z__[i__ - 1] = 1. / sl;
    if (z__[i__ - 1] <= scbd)
    {
      goto L189;
    }
    sl = 1. / scbd;
    z__[i__ - 1] = scbd;
    L189: i__1 =  *n;
    for (j = 1; j <= i__1; ++j)
    {
      /* L190: */
      q_1.v[i__ + j * 20-21] = sl * q_1.v[i__ + j * 20-21];
    }
    /* L195: */
  }

  /* .....CALCULATE A NEW SET OF ORTHOGONAL DIRECTIONS BEFORE REPEATING */
  /*     THE MAIN LOOP. */
  /*     FIRST TRANSPOSE V FOR MINFIT: */

  L200: i__2 =  *n;
  for (i__ = 2; i__ <= i__2; ++i__)
  {
    im1 = i__ - 1;
    i__1 = im1;
    for (j = 1; j <= i__1; ++j)
    {
      s = q_1.v[i__ + j * 20-21];
      q_1.v[i__ + j * 20-21] = q_1.v[j + i__ * 20-21];
      /* L210: */
      q_1.v[j + i__ * 20-21] = s;
    }
    /* L220: */
  }

  /* .....CALL MINFIT TO FIND THE SINGULAR VALUE DECOMPOSITION OF V. */
  /*     THIS GIVES THE PRINCIPAL VALUES AND PRINCIPAL DIRECTIONS OF THE */
  /*     APPROXIMATING QUADRATIC FORM WITHOUT SQUARING THE CONDITION */
  /*     NUMBER..... */

  minfit_(&idim, n, machep, &vsmall, q_1.v, d__);

  /* .....UNSCALE THE AXES..... */

  if (scbd <= 1.)
  {
    goto L250;
  }
  i__2 =  *n;
  for (i__ = 1; i__ <= i__2; ++i__)
  {
    s = z__[i__ - 1];
    i__1 =  *n;
    for (j = 1; j <= i__1; ++j)
    {
      /* L225: */
      q_1.v[i__ + j * 20-21] = s * q_1.v[i__ + j * 20-21];
    }
    /* L230: */
  }
  i__2 =  *n;
  for (i__ = 1; i__ <= i__2; ++i__)
  {
    s = 0.;
    i__1 =  *n;
    for (j = 1; j <= i__1; ++j)
    {
      /* L235: */
      /* Computing 2nd power */
      d__1 = q_1.v[j + i__ * 20-21];
      s += d__1 * d__1;
    }
    s = System::Math::Sqrt(s);
    d__[i__ - 1] = s * d__[i__ - 1];
    s = 1 / s;
    i__1 =  *n;
    for (j = 1; j <= i__1; ++j)
    {
      /* L240: */
      q_1.v[j + i__ * 20-21] = s * q_1.v[j + i__ * 20-21];
    }
    /* L245: */
  }

  L250: i__2 =  *n;
  for (i__ = 1; i__ <= i__2; ++i__)
  {
    dni = dn * d__[i__ - 1];
    if (dni > large)
    {
      goto L265;
    }
    if (dni < small)
    {
      goto L260;
    }
    d__[i__ - 1] = 1 / (dni *dni);
    goto L270;
    L260: d__[i__ - 1] = vlarge;
    goto L270;
    L265: d__[i__ - 1] = vsmall;
    L270: ;
  }

  /* .....SORT THE EIGENVALUES AND EIGENVECTORS..... */

  sort_(&idim, n, d__, q_1.v);
  global_1.dmin__ = d__[ *n - 1];
  if (global_1.dmin__ < small)
  {
    global_1.dmin__ = small;
  }
  illc = FALSE_;
  if (m2 *d__[0] > global_1.dmin__)
  {
    illc = TRUE_;
  }
  /* if (*prin > 1 && scbd > 1.) {
  vcprnt_(&c__2, z__, n);
  }
  if (*prin > 1) {
  vcprnt_(&c__3, d__, n);
  }
  if (*prin > 3) {
  maprnt_(&c__2, q_1.v, &idim, n);
  }*/
  /* .....THE MAIN LOOP ENDS HERE..... */

  goto L40;

  /* .....RETURN..... */

  L400: 
  /*if (*prin > 0) {
  vcprnt_(&c__4, &x[1], n);
  }*/
  ret_val = global_1.fx;
  return ret_val;
} /* praxis_ */

 /* Subroutine */int SmartMathLibrary::InternalRoutines
   ::PrincipalAxisOptimizerProvider::minfit_(integer *m, integer *n, doublereal
   *machep, doublereal *tol, doublereal *ab, doublereal *q)
{
  /* Format strings */
  static char fmt_1000[] = "(\002 QR FAILED\002)";

  /* System generated locals */
  integer ab_dim1, ab_offset, i__1, i__2, i__3;
  doublereal d__1, d__2;

  /* Local variables */
  static doublereal c__, e[20], f, g, h__;
  static integer i__, j, k, l;
  static doublereal s, x, y, z__;
  static integer l2, ii, kk, kt, ll2, lp1;
  static doublereal eps, temp;

  /* Fortran I/O blocks */
  static cilist io___52 = 
  {
    0, 6, 0, fmt_1000, 0
  };


  /* ...AN IMPROVED VERSION OF MINFIT (SEE GOLUB AND REINSCH, 1969) */
  /*   RESTRICTED TO M=N,P=0. */
  /*   THE SINGULAR VALUES OF THE ARRAY AB ARE RETURNED IN Q AND AB IS */
  /*   OVERWRITTEN WITH THE ORTHOGONAL MATRIX V SUCH THAT U.DIAG(Q) = AB.V, */
  /*   WHERE U IS ANOTHER ORTHOGONAL MATRIX. */
  /* ...HOUSEHOLDER'S REDUCTION TO BIDIAGONAL FORM... */
  /* Parameter adjustments */
  --q;
  ab_dim1 =  *m;
  ab_offset = 1+ab_dim1;
  ab -= ab_offset;

  /* Function Body */
  if (*n == 1)
  {
    goto L200;
  }
  eps =  *machep;
  g = 0.;
  x = 0.;
  i__1 =  *n;
  for (i__ = 1; i__ <= i__1; ++i__)
  {
    e[i__ - 1] = g;
    s = 0.;
    l = i__ + 1;
    i__2 =  *n;
    for (j = i__; j <= i__2; ++j)
    {
      /* L1: */
      /* Computing 2nd power */
      d__1 = ab[j + i__ * ab_dim1];
      s += d__1 * d__1;
    }
    g = 0.;
    if (s <  *tol)
    {
      goto L4;
    }
    f = ab[i__ + i__ * ab_dim1];
    g = System::Math::Sqrt(s);
    if (f >= 0.)
    {
      g =  - g;
    }
    h__ = f * g - s;
    ab[i__ + i__ * ab_dim1] = f - g;
    if (l >  *n)
    {
      goto L4;
    }
    i__2 =  *n;
    for (j = l; j <= i__2; ++j)
    {
      f = 0.;
      i__3 =  *n;
      for (k = i__; k <= i__3; ++k)
      {
        /* L2: */
        f += ab[k + i__ * ab_dim1] *ab[k + j * ab_dim1];
      }
      f /= h__;
      i__3 =  *n;
      for (k = i__; k <= i__3; ++k)
      {
        /* L3: */
        ab[k + j * ab_dim1] += f * ab[k + i__ * ab_dim1];
      }
    }
    L4: q[i__] = g;
    s = 0.;
    if (i__ ==  *n)
    {
      goto L6;
    }
    i__3 =  *n;
    for (j = l; j <= i__3; ++j)
    {
      /* L5: */
      s += ab[i__ + j * ab_dim1] *ab[i__ + j * ab_dim1];
    }
    L6: g = 0.;
    if (s <  *tol)
    {
      goto L10;
    }
    if (i__ ==  *n)
    {
      goto L16;
    }
    f = ab[i__ + (i__ + 1) *ab_dim1];
    L16: g = System::Math::Sqrt(s);
    if (f >= 0.)
    {
      g =  - g;
    }
    h__ = f * g - s;
    if (i__ ==  *n)
    {
      goto L10;
    }
    ab[i__ + (i__ + 1) *ab_dim1] = f - g;
    i__3 =  *n;
    for (j = l; j <= i__3; ++j)
    {
      /* L7: */
      e[j - 1] = ab[i__ + j * ab_dim1] / h__;
    }
    i__3 =  *n;
    for (j = l; j <= i__3; ++j)
    {
      s = 0.;
      i__2 =  *n;
      for (k = l; k <= i__2; ++k)
      {
        /* L8: */
        s += ab[j + k * ab_dim1] *ab[i__ + k * ab_dim1];
      }
      i__2 =  *n;
      for (k = l; k <= i__2; ++k)
      {
        /* L9: */
        ab[j + k * ab_dim1] += s * e[k - 1];
      }
    }
    L10: y = (d__1 = q[i__], abs(d__1)) + (d__2 = e[i__ - 1], abs(d__2));
    /* L11: */
    if (y > x)
    {
      x = y;
    }
  }
  /* ...ACCUMULATION OF RIGHT-HAND TRANSFORMATIONS... */
  ab[ *n +  *n * ab_dim1] = 1.;
  g = e[ *n - 1];
  l =  *n;
  i__1 =  *n;
  for (ii = 2; ii <= i__1; ++ii)
  {
    i__ =  *n - ii + 1;
    if (g == 0.)
    {
      goto L23;
    }
    h__ = ab[i__ + (i__ + 1) *ab_dim1] *g;
    i__2 =  *n;
    for (j = l; j <= i__2; ++j)
    {
      /* L20: */
      ab[j + i__ * ab_dim1] = ab[i__ + j * ab_dim1] / h__;
    }
    i__2 =  *n;
    for (j = l; j <= i__2; ++j)
    {
      s = 0.;
      i__3 =  *n;
      for (k = l; k <= i__3; ++k)
      {
        /* L21: */
        s += ab[i__ + k * ab_dim1] *ab[k + j * ab_dim1];
      }
      i__3 =  *n;
      for (k = l; k <= i__3; ++k)
      {
        /* L22: */
        ab[k + j * ab_dim1] += s * ab[k + i__ * ab_dim1];
      }
    }
    L23: i__3 =  *n;
    for (j = l; j <= i__3; ++j)
    {
      ab[i__ + j * ab_dim1] = 0.;
      /* L24: */
      ab[j + i__ * ab_dim1] = 0.;
    }
    ab[i__ + i__ * ab_dim1] = 1.;
    g = e[i__ - 1];
    /* L25: */
    l = i__;
  }
  /* ...DIAGONALIZATION OF THE BIDIAGONAL FORM... */
  /* L100: */
  eps *= x;
  i__1 =  *n;
  for (kk = 1; kk <= i__1; ++kk)
  {
    k =  *n - kk + 1;
    kt = 0;
    L101: ++kt;
    if (kt <= 30)
    {
      goto L102;
    }
    e[k - 1] = 0.;
    //s_wsfe(&io___52);
    //e_wsfe();
    L102: i__3 = k;
    for (ll2 = 1; ll2 <= i__3; ++ll2)
    {
      l2 = k - ll2 + 1;
      l = l2;
      if ((d__1 = e[l - 1], abs(d__1)) <= eps)
      {
        goto L120;
      }
      if (l == 1)
      {
        goto L103;
      }
      if ((d__1 = q[l - 1], abs(d__1)) <= eps)
      {
        goto L110;
      }
      L103: ;
    }
    /* ...CANCELLATION OF E(L) IF L>1... */
    L110: c__ = 0.;
    s = 1.;
    i__3 = k;
    for (i__ = l; i__ <= i__3; ++i__)
    {
      f = s * e[i__ - 1];
      e[i__ - 1] = c__ * e[i__ - 1];
      if (abs(f) <= eps)
      {
        goto L120;
      }
      g = q[i__];
      /* ...Q(I) = H = DSystem::Math::Sqrt(G*G + F*F)... */
      if (abs(f) < abs(g))
      {
        goto L113;
      }
      if (f != 0.)
      {
        goto L112;
      }
      else
      {
        goto L111;
      }
      L111: h__ = 0.;
      goto L114;
      L112: 
      /* Computing 2nd power */
      d__1 = g / f;
      h__ = abs(f) *System::Math::Sqrt(d__1 *d__1 + 1);
      goto L114;
      L113: 
      /* Computing 2nd power */
      d__1 = f / g;
      h__ = abs(g) *System::Math::Sqrt(d__1 *d__1 + 1);
      L114: q[i__] = h__;
      if (h__ != 0.)
      {
        goto L115;
      }
      g = 1.;
      h__ = 1.;
      L115: c__ = g / h__;
      /* L116: */
      s =  - f / h__;
    }
    /* ...TEST FOR CONVERGENCE... */
    L120: z__ = q[k];
    if (l == k)
    {
      goto L140;
    }
    /* ...SHIFT FROM BOTTOM 2*2 MINOR... */
    x = q[l];
    y = q[k - 1];
    g = e[k - 2];
    h__ = e[k - 1];
    f = ((y - z__)*(y + z__) + (g - h__)*(g + h__)) / (h__ *2 * y);
    g = System::Math::Sqrt(f *f + 1.);
    temp = f - g;
    if (f >= 0.)
    {
      temp = f + g;
    }
    f = ((x - z__)*(x + z__) + h__ *(y / temp - h__)) / x;
    /* ...NEXT QR TRANSFORMATION... */
    c__ = 1.;
    s = 1.;
    lp1 = l + 1;
    if (lp1 > k)
    {
      goto L133;
    }
    i__3 = k;
    for (i__ = lp1; i__ <= i__3; ++i__)
    {
      g = e[i__ - 1];
      y = q[i__];
      h__ = s * g;
      g *= c__;
      if (abs(f) < abs(h__))
      {
        goto L123;
      }
      if (f != 0.)
      {
        goto L122;
      }
      else
      {
        goto L121;
      }
      L121: z__ = 0.;
      goto L124;
      L122: 
      /* Computing 2nd power */
      d__1 = h__ / f;
      z__ = abs(f) *System::Math::Sqrt(d__1 *d__1 + 1);
      goto L124;
      L123: 
      /* Computing 2nd power */
      d__1 = f / h__;
      z__ = abs(h__) *System::Math::Sqrt(d__1 *d__1 + 1);
      L124: e[i__ - 2] = z__;
      if (z__ != 0.)
      {
        goto L125;
      }
      f = 1.;
      z__ = 1.;
      L125: c__ = f / z__;
      s = h__ / z__;
      f = x * c__ + g * s;
      g =  - x * s + g * c__;
      h__ = y * s;
      y *= c__;
      i__2 =  *n;
      for (j = 1; j <= i__2; ++j)
      {
        x = ab[j + (i__ - 1) *ab_dim1];
        z__ = ab[j + i__ * ab_dim1];
        ab[j + (i__ - 1) *ab_dim1] = x * c__ + z__ * s;
        /* L126: */
        ab[j + i__ * ab_dim1] =  - x * s + z__ * c__;
      }
      if (abs(f) < abs(h__))
      {
        goto L129;
      }
      if (f != 0.)
      {
        goto L128;
      }
      else
      {
        goto L127;
      }
      L127: z__ = 0.;
      goto L130;
      L128: 
      /* Computing 2nd power */
      d__1 = h__ / f;
      z__ = abs(f) *System::Math::Sqrt(d__1 *d__1 + 1);
      goto L130;
      L129: 
      /* Computing 2nd power */
      d__1 = f / h__;
      z__ = abs(h__) *System::Math::Sqrt(d__1 *d__1 + 1);
      L130: q[i__ - 1] = z__;
      if (z__ != 0.)
      {
        goto L131;
      }
      f = 1.;
      z__ = 1.;
      L131: c__ = f / z__;
      s = h__ / z__;
      f = c__ * g + s * y;
      /* L132: */
      x =  - s * g + c__ * y;
    }
    L133: e[l - 1] = 0.;
    e[k - 1] = f;
    q[k] = x;
    goto L101;
    /* ...CONVERGENCE:  Q(K) IS MADE NON-NEGATIVE... */
    L140: if (z__ >= 0.)
    {
      goto L150;
    }
    q[k] =  - z__;
    i__3 =  *n;
    for (j = 1; j <= i__3; ++j)
    {
      /* L141: */
      ab[j + k * ab_dim1] =  - ab[j + k * ab_dim1];
    }
    L150: ;
  }
  return 0;
  L200: q[1] = ab[ab_dim1 + 1];
  ab[ab_dim1 + 1] = 1.;
  return 0;
} /* minfit_ */

 /* Subroutine */int SmartMathLibrary::InternalRoutines
   ::PrincipalAxisOptimizerProvider::min_(integer *n, integer *j, integer *nits,
   doublereal *d2, doublereal *x1, doublereal *f1, logical *fk,
   PrincipalAxisHelperFunction ^ f, doublereal *x, doublereal *t, doublereal
   *machep, doublereal *h__)
{
  /* System generated locals */
  integer i__1;
  doublereal d__1, d__2;

  /* Local variables */
  static integer i__, k;
  static doublereal s, d1, f0, f2, m2, m4, t2, x2, fm;
  static logical dz;
  static doublereal xm, sf1, sx1;
  static doublereal temp, small;

  /* ...THE SUBROUTINE MIN MINIMIZES F FROM X IN THE DIRECTION V(*,J) UNLESS */
  /*   J IS LESS THAN 1, WHEN A QUADRATIC SEARCH IS MADE IN THE PLANE */
  /*   DEFINED BY Q0,Q1,X. */
  /*   D2 IS EITHER ZERO OR AN APPROXIMATION TO HALF F". */
  /*   ON ENTRY, X1 IS AN ESTIMATE OF THE DISTANCE FROM X TO THE MINIMUM */
  /*   ALONG V(*,J) (OR, IF J=0, A CURVE).  ON RETURN, X1 IS THE DISTANCE */
  /*   FOUND. */
  /*   IF FK=.TRUE., THEN F1 IS FLIN(X1).  OTHERWISE X1 AND F1 ARE IGNORED */
  /*   ON ENTRY UNLESS FINAL FX IS GREATER THAN F1. */
  /*   NITS CONTROLS THE NUMBER OF TIMES AN ATTEMPT WILL BE MADE TO HALVE */
  /*   THE INTERVAL. */
  /* Parameter adjustments */
  --x;

  /* Function Body */
  /* Computing 2nd power */
  d__1 =  *machep;
  small = d__1 * d__1;
  m2 = System::Math::Sqrt(*machep);
  m4 = System::Math::Sqrt(m2);
  sf1 =  *f1;
  sx1 =  *x1;
  k = 0;
  xm = 0.;
  fm = global_1.fx;
  f0 = global_1.fx;
  dz =  *d2 <  *machep;
  /* ...FIND THE STEP SIZE... */
  s = 0.;
  i__1 =  *n;
  for (i__ = 1; i__ <= i__1; ++i__)
  {
    /* L1: */
    /* Computing 2nd power */
    d__1 = x[i__];
    s += d__1 * d__1;
  }
  s = System::Math::Sqrt(s);
  temp =  *d2;
  if (dz)
  {
    temp = global_1.dmin__;
  }
  t2 = m4 * System::Math::Sqrt(abs(global_1.fx) / temp + s * global_1.ldt) + m2
    * global_1.ldt;
  s = m4 * s +  *t;
  if (dz && t2 > s)
  {
    t2 = s;
  }
  t2 = max(t2, small);
  /* Computing MIN */
  d__1 = t2, d__2 =  *h__ * .01;
  t2 = min(d__1, d__2);
  if (!(*fk) ||  *f1 > fm)
  {
    goto L2;
  }
  xm =  *x1;
  fm =  *f1;
  L2: if (*fk && abs(*x1) >= t2)
  {
    goto L3;
  }
  temp = 1.;
  if (*x1 < 0.)
  {
    temp =  - 1.;
  }
  *x1 = temp * t2;
  *f1 = flin_(n, j, x1, f, &x[1], &global_1.nf);
  L3: if (*f1 > fm)
  {
    goto L4;
  }
  xm =  *x1;
  fm =  *f1;
  L4: if (!dz)
  {
    goto L6;
  }
  /* ...EVALUATE FLIN AT ANOTHER POINT AND ESTIMATE THE SECOND DERIVATIVE... */
  x2 =  - (*x1);
  if (f0 >=  *f1)
  {
    x2 =  *x1 * 2.;
  }
  f2 = flin_(n, j, &x2, f, &x[1], &global_1.nf);
  if (f2 > fm)
  {
    goto L5;
  }
  xm = x2;
  fm = f2;
  L5: *d2 = (x2 *(*f1 - f0) -  *x1 *(f2 - f0)) / (*x1 * x2 *(*x1 - x2));
  /* ...ESTIMATE THE FIRST DERIVATIVE AT 0... */
  L6: d1 = (*f1 - f0) /  *x1 -  *x1 **d2;
  dz = TRUE_;
  /* ...PREDICT THE MINIMUM... */
  if (*d2 > small)
  {
    goto L7;
  }
  x2 =  *h__;
  if (d1 >= 0.)
  {
    x2 =  - x2;
  }
  goto L8;
  L7: x2 = d1 *  - .5 /  *d2;
  L8: if (abs(x2) <=  *h__)
  {
    goto L11;
  }
  if (x2 <= 0.)
  {
    goto L9;
  }
  else
  {
    goto L10;
  }
  L9: x2 =  - (*h__);
  goto L11;
  L10: x2 =  *h__;
  /* ...EVALUATE F AT THE PREDICTED MINIMUM... */
  L11: f2 = flin_(n, j, &x2, f, &x[1], &global_1.nf);
  if (k >=  *nits || f2 <= f0)
  {
    goto L12;
  }
  /* ...NO SUCCESS, SO TRY AGAIN... */
  ++k;
  if (f0 <  *f1 &&  *x1 * x2 > 0.)
  {
    goto L4;
  }
  x2 *= .5;
  goto L11;
  /* ...INCREMENT THE ONE-DIMENSIONAL SEARCH COUNTER... */
  L12: ++global_1.nl;
  if (f2 <= fm)
  {
    goto L13;
  }
  x2 = xm;
  goto L14;
  L13: fm = f2;
  /* ...GET A NEW ESTIMATE OF THE SECOND DERIVATIVE... */
  L14: if ((d__1 = x2 *(x2 -  *x1), abs(d__1)) <= small)
  {
    goto L15;
  }
  *d2 = (x2 *(*f1 - f0) -  *x1 *(fm - f0)) / (*x1 * x2 *(*x1 - x2));
  goto L16;
  L15: if (k > 0)
  {
    *d2 = 0.;
  }
  L16: if (*d2 <= small)
  {
    *d2 = small;
  }
  *x1 = x2;
  global_1.fx = fm;
  if (sf1 >= global_1.fx)
  {
    goto L17;
  }
  global_1.fx = sf1;
  *x1 = sx1;
  /* ...UPDATE X FOR LINEAR BUT NOT PARABOLIC SEARCH... */
  L17: if (*j == 0)
  {
    return 0;
  }
  i__1 =  *n;
  for (i__ = 1; i__ <= i__1; ++i__)
  {
    /* L18: */
    x[i__] +=  *x1 * q_1.v[i__ +  *j * 20-21];
  }
  return 0;
} /* min_ */

doublereal SmartMathLibrary::InternalRoutines::PrincipalAxisOptimizerProvider
  ::flin_(integer *n, integer *j, doublereal *l, PrincipalAxisHelperFunction ^
  f, doublereal *x, integer *nf)
{
  /* System generated locals */
  integer i__1;
  doublereal ret_val;

  /* Local variables */
  static integer i__;
  static doublereal t[20];

  /* ...FLIN IS THE FUNCTION OF ONE REAL VARIABLE L THAT IS MINIMIZED */
  /*   BY THE SUBROUTINE MIN... */
  /* Parameter adjustments */
  --x;

  /* Function Body */
  if (*j == 0)
  {
    goto L2;
  }
  /* ...THE SEARCH IS LINEAR... */
  i__1 =  *n;
  for (i__ = 1; i__ <= i__1; ++i__)
  {
    /* L1: */
    t[i__ - 1] = x[i__] +  *l * q_1.v[i__ +  *j * 20-21];
  }
  goto L4;
  /* ...THE SEARCH IS ALONG A PARABOLIC SPACE CURVE... */
  L2: q_1.qa =  *l *(*l - q_1.qd1) / (q_1.qd0 *(q_1.qd0 + q_1.qd1));
  q_1.qb = (*l + q_1.qd0)*(q_1.qd1 -  *l) / (q_1.qd0 *q_1.qd1);
  q_1.qc =  *l *(*l + q_1.qd0) / (q_1.qd1 *(q_1.qd0 + q_1.qd1));
  i__1 =  *n;
  for (i__ = 1; i__ <= i__1; ++i__)
  {
    /* L3: */
    t[i__ - 1] = q_1.qa * q_1.q0[i__ - 1] + q_1.qb * x[i__] + q_1.qc *
      q_1.q1[i__ - 1];
  }
  /* ...THE FUNCTION EVALUATION COUNTER NF IS INCREMENTED... */
  L4: ++(*nf);
  ret_val = f(t, n);
  return ret_val;
} /* flin_ */

 /* Subroutine */int SmartMathLibrary::InternalRoutines
   ::PrincipalAxisOptimizerProvider::sort_(integer *m, integer *n, doublereal
   *d__, doublereal *v)
{
  /* System generated locals */
  integer v_dim1, v_offset, i__1, i__2;

  /* Local variables */
  static integer i__, j, k;
  static doublereal s;
  static integer ip1, nm1;

  /* ...SORTS THE ELEMENTS OF D(N) INTO DESCENDING ORDER AND MOVES THE */
  /*   CORRESPONDING COLUMNS OF V(N,N). */
  /*   M IS THE ROW DIMENSION OF V AS DECLARED IN THE CALLING PROGRAM. */
  /* Parameter adjustments */
  v_dim1 =  *m;
  v_offset = 1+v_dim1;
  v -= v_offset;
  --d__;

  /* Function Body */
  if (*n == 1)
  {
    return 0;
  }
  nm1 =  *n - 1;
  i__1 = nm1;
  for (i__ = 1; i__ <= i__1; ++i__)
  {
    k = i__;
    s = d__[i__];
    ip1 = i__ + 1;
    i__2 =  *n;
    for (j = ip1; j <= i__2; ++j)
    {
      if (d__[j] <= s)
      {
        goto L1;
      }
      k = j;
      s = d__[j];
      L1: ;
    }
    if (k <= i__)
    {
      goto L3;
    }
    d__[k] = d__[i__];
    d__[i__] = s;
    i__2 =  *n;
    for (j = 1; j <= i__2; ++j)
    {
      s = v[j + i__ * v_dim1];
      v[j + i__ * v_dim1] = v[j + k * v_dim1];
      /* L2: */
      v[j + k * v_dim1] = s;
    }
    L3: ;
  }
  return 0;
} /* sort_ */

 /* Subroutine */int SmartMathLibrary::InternalRoutines
   ::PrincipalAxisOptimizerProvider::quad_(integer *n,
   PrincipalAxisHelperFunction ^ f, doublereal *x, doublereal *t, doublereal
   *machep, doublereal *h__)
{
  /* System generated locals */
  integer i__1;
  doublereal d__1;

  /* Local variables */
  static integer i__;
  static doublereal l, s;
  static doublereal value;

  /* ...QUAD LOOKS FOR THE MINIMUM OF F ALONG A CURVE DEFINED BY Q0,Q1,X... */
  /* Parameter adjustments */
  --x;

  /* Function Body */
  s = global_1.fx;
  global_1.fx = q_1.qf1;
  q_1.qf1 = s;
  q_1.qd1 = 0.;
  i__1 =  *n;
  for (i__ = 1; i__ <= i__1; ++i__)
  {
    s = x[i__];
    l = q_1.q1[i__ - 1];
    x[i__] = l;
    q_1.q1[i__ - 1] = s;
    /* L1: */
    /* Computing 2nd power */
    d__1 = s - l;
    q_1.qd1 += d__1 * d__1;
  }
  q_1.qd1 = System::Math::Sqrt(q_1.qd1);
  l = q_1.qd1;
  s = 0.;
  if (q_1.qd0 <= 0. || q_1.qd1 <= 0. || global_1.nl <  *n * 3 **n)
  {
    goto L2;
  }
  value = q_1.qf1;
  min_(n, &c__0, &c__2, &s, &l, &value, &c_true, f, &x[1], t, machep, h__);
  q_1.qa = l *(l - q_1.qd1) / (q_1.qd0 *(q_1.qd0 + q_1.qd1));
  q_1.qb = (l + q_1.qd0)*(q_1.qd1 - l) / (q_1.qd0 *q_1.qd1);
  q_1.qc = l *(l + q_1.qd0) / (q_1.qd1 *(q_1.qd0 + q_1.qd1));
  goto L3;
  L2: global_1.fx = q_1.qf1;
  q_1.qa = 0.;
  q_1.qb = q_1.qa;
  q_1.qc = 1.;
  L3: q_1.qd0 = q_1.qd1;
  i__1 =  *n;
  for (i__ = 1; i__ <= i__1; ++i__)
  {
    s = q_1.q0[i__ - 1];
    q_1.q0[i__ - 1] = x[i__];
    /* L4: */
    x[i__] = q_1.qa * s + q_1.qb * x[i__] + q_1.qc * q_1.q1[i__ - 1];
  }
  return 0;
} /* quad_ */

doublereal SmartMathLibrary::InternalRoutines::PrincipalAxisOptimizerProvider
  ::random_(integer *naught)
{
  /* Initialized data */

  static logical init = FALSE_;

  /* System generated locals */
  doublereal ret_val;

  /* Local variables */
  static integer i__, j, q, r__;
  static doublereal ran1;
  static integer ran2;
  static doublereal ran3[127], half;

  if (init)
  {
    goto L3;
  }
  r__ =  *naught % 8190+1;
  ran2 = 128;
  for (i__ = 1; i__ <= 127; ++i__)
  {
    --ran2;
    ran1 =  - 36028797018963968.;
    for (j = 1; j <= 7; ++j)
    {
      r__ = r__ * 1756 % 8191;
      q = r__ / 32;
      /* L1: */
      ran1 = (ran1 + q) *.00390625;
    }
    /* L2: */
    ran3[ran2 - 1] = ran1;
  }
  init = TRUE_;
  L3: if (ran2 == 1)
  {
    ran2 = 128;
  }
  --ran2;
  ran1 += ran3[ran2 - 1];
  half = .5;
  if (ran1 >= 0.)
  {
    half =  - half;
  }
  ran1 += half;
  ran3[ran2 - 1] = ran1;
  ret_val = ran1 + .5;
  return ret_val;
} /* random_ */
integer SmartMathLibrary::InternalRoutines::PrincipalAxisOptimizerProvider::pow_ii
  (integer *ap, integer *bp)
{
  integer pow, x, n;
  unsigned long u;

  x =  *ap;
  n =  *bp;

  if (n <= 0)
  {
    if (n == 0 || x == 1)
      return 1;
    if (x !=  - 1)
      return x == 0 ? 1 / x: 0;
    n =  - n;
  }
  u = n;
  for (pow = 1;;)
  {
    if (u &01)
      pow *= x;
    if (u >>= 1)
      x *= x;
    else
      break;
  }
  return (pow);
}

double SmartMathLibrary::InternalRoutines::PrincipalAxisOptimizerProvider::d_lg10
  (doublereal *x)
{
  return (0.43429448190325182765 *System::Math::Log(*x));
}
